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Abstract 



In this paper we present in one-dimensional space a numerical solution of a partial differential equation of fractional order. This 
equation describes a process of anomalous diffusion. The process arises from the interactions within the complex and non-homogeneous 
background. We presented a numerical method which bases on the finite differences method. We considered pure initial and boundary- 
initial value problems for the equation with the Riesz-Feller fractional derivative. In the final part of this paper sample results of 
simulation were shown. 



Keywords: anomalous diffusion, fractional calculus, Riesz-Feller derivative, finite difference method, boundary value problem 



1. Introduction 

Anomalous diffusion is a phenomenon strongly connected 
with the interactions within complex and non-homogeneous back- 
ground. This phenomenon is observed in transport of fluid in 
porous materials, in the chaotic heat baths, amorphous semicon- 
ductors, particle dynamics inside polymer network, two-dimen- 
sional rotating flow and also in econophysics. Phenomenon of 
anomalous diffusion deviates from the standard diffusion behaviour. 
In opposite to standard diffusion where linear form in the mean 
square displacement ( x 2 (t)) ~ kit of diffusing particle over 
time occurs, anomalous diffusion is characterized by the non- 
linear one ( x 2 (t)) ~ fc 7 i 7 , for 7 6 (0, 2]. In this phenomenon 
may exist dependence ( x 2 (t)) — > 00 , which is characterized by 
occurrence of rare but extremely large jumps of diffusing particle 
- well-known as the Levy motion or the Levy flights. Ordinary 
diffusion follows Gaussian statistics and Fick's second law for 
finding running process at time t whereas anomalous diffusion 
follows non-Gaussian statistic or can be interpreted as the Levy 
stable densities. 

Many authors proposed models which base on linear and non- 
linear forms of differential equations. Such models can simulate 
anomalous diffusion but they don't reflect its real behaviour. Sev- 
eral authors |2| 151171 151 II II Il3l 1151 apply fractional calculus in 
modelling of this type of diffusion. This means that time and 
spatial derivatives in the classical diffusion equation are replaced 
by fractional ones. In comparison to derivatives of integer order, 
which depend on the local behaviour of the function, derivatives 
of fractional order accumulate the whole history of this function. 

2. Mathematical background 



unit of measure [m a /s]. According to (7| ll II the Riesz-Feller 
fractional operator for 0<a<2, a^l for one-variable func- 
tion u{x) is 

——^u(x) = x Dgu (x) = - [c L (a, 6) - x D"u (x) 

a\x\ (2) 



where 
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dx 



+c R {a, 9) xD+^u (x)] . 
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and coefficients cl (ct, 9), cr (a, 9) (for < a < 2, q / 1, and 
for \9\ < min (a, 2 — a)), are defined as 



(5) 



The fractional operators of order a: - x I"u (x) and xl^u (x) 
are defined as the left- and right-side of Weyl fractional integrals 
|6|I71 |14I|15II16I which definitions are 



J"u(x) = 



«(0 



r(o) J_ x (x - 



— d£, 

1 — ot ^' 



(6) 



In this paper, we consider an equation in the following form 

d d a 

— C(x,t) — k a ^-r—r^C(x,t), !>0,l6l, (1) 

at d\x\ 

where C(x, t) is a field variable, j$^m C(x, t) is the Riesz-Feller 
fractional operator II 311 161 . a is the real order of this operator, k a 
is the coefficient of generalized (anomalous) diffusion with the 



Considering Eqn (1) we obtain the classical diffusion equa- 
tion for a = 2, i.e. the heat transfer equation. If a — 1, and the 
parameter of skewness 9 admits extreme values in ^5}, the trans- 
port equation is noted. Therefore we assume variations of the 
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parameter a within the range < a < 2. Analysing behaviour 
of the parameter a < 2 in Eqn (1), we found some combination 
between transport and propagation processes. 

For analytic solution of Eqn (1) we can apply Green func- 
tions |6|. We numerically solve Eqn (1) when additional non- 
linear term may occur. Some numerical methods used in solution 
of fractional differential equations can be found in |7 1. However 
they apply the infinite domain. 

In this work we will consider Eqn (1) limited for 1 < a < 2 
in one dimensional domain fl : L < x < R with the boundary- 
value conditions of the first kind (the Dirichlet conditions) as 



L 



C(L,t)=g L (t). 



t > 0, 



(8) 



(9) 



[ x = R: C (R,t) = g R (t) , 
and with the initial-value condition 

C(x,t)\ t=0 = c o («) ■ 
3. Numerical method 



According to the finite difference method fTll4*ll5ll9l ll2l we 
consider a discrete from of Eqn (1) both in time and space. In the 
previous work 1 3 1 we solved numerically the anomalous diffusion 
equation similar to the Eqn ( 1 ) with the time-fractional derivative. 
We called this method FFDM (Fractional FDM). The problem of 
solving of Eqn (1) lies in properly approximation of the Riesz- 
Feller derivative (2) in numerical scheme. 

3. 1. Approximation of the Riesz-Feller derivative 

We begin numerical analysis from discrete forms of opera- 
tors (6) and (7). We introduce homogenous spatial grid — oo < 
. . . < Xi-2 < Xi-l < Xi < Xi+l < X i+ 2 < . . . < oo with 
the step h = x k — Xk-i and we denote value of function u in 
the point Xk as u& = u (xt), for fc £ Z. In order to simplify no- 
tations we take here the function of one variable. For numerical 
integration scheme we assumed the trapezoidal rule. The inte- 
gral (6) in point Xi of the grid is replaced by the sum of discrete 
integrals as 



r(a 



1 °° 
W) E 



■(£) 



1 — a ^ 



(10) 



and using linear interpolation of function u in every sub-interval 

[Xi-k-l,Xi-k] 
* fr\ u i-k — Ui-k-lXi-k — Ui-kXi-k-1 

u (0 = ~ h C + ~ h (ID 

we have 

1 °° 



«* (0 



r(a)^ J (*i-0 
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— d ? 



^ oo 
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where 
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(13) 
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After next transforms we can write 

oo 

where 



(a) 



(a) 

= 



X 



fe r(2 + a) 
1, 

1+a 2k 1+a 



(15) 



(16) 



(fc + l) J 



for k = 0, 
+ (fc-l) 1+Q ,forfc = l,...,oo. 



Similar to previous considerations we approximate operator 
i/^u (x) in the point x% and finally we obtain 



rc* _ i a \ ^ (a) 



(17) 



where coefficients have identical forms as (16). 

In the next step we analyse operator (2). It can be expressed 
in the form (in order to simplify this we denote cl = cl (a, 9) 

and cr = c R (a, 9) ) 



c Dqu (x) = 
d 2 r 



(18) 



CL 



dx 2 



Jl a U (x)] + Cr — [ X I+J?U(. 



We used the central difference scheme for the second spatial deriva- 
tive in the point x, and we obtain 



(19) 



72 fl r\ rl — OL I r2 — CK 
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/i 2 

r2 — a„ 
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After numerous transforms we obtain the final form as 



i E u ^ 



(a) 



(20) 



where coefficients are 



(a) 
W k = 



-1 



r(4-a) 

• ((|fc| + 2) 3 " Q - 4 (|fc| + I) 3 "" + 6 \kf~ a 

-4 (\k\ - lf- a + (\k\ - 2) 3 ~ Q ) cl, for k < -2 
( 3 3-« _2 5 " Q +6) c L + cr, for ft = -1 

(2 3 " a - 4) (cl + c«) , for k = . 

( 3 s-« „ 2 b ' a +6) c fl + c L , for ft = 1 

((k + 2°f- a -4(k+ l) 3 ~ a + 6k 3 ~ a 

~4(k~l)' A ' a + {k~2f~ a )c R , forfc>2 

Assuming a = 2 and 6 = we have cl (2, 0) = cr (2, 0) 
- 1 and we obtain 



(21) 



(2) 



f o, 


for k < -2 


1, 


forfc = -1 


-2, 


for fc = 


1, 


for ft = 1 


, °. 


for fc > 2 



(22) 



These coeeficients are identical as for wide known the central 
difference scheme for the second derivative. Also when a — > 1 + 
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and 8 — after arduous calculations of limits we obtain coeffi- 
cients 



(i+) 1 
wl = — X 



hi 



2tt 

(|fc| + l) 4(|fc|+1)2 



(23) 



x4(|fc|-ir 



(|fc|+2) (|fc|+2)2 |fc| 6fc2 (|fc|-2) (|fc| - 2)2 
16 In 2 - 9 In 3, 
-8 In 2, 

16 In 2 - 9 In 3, 

(fc+l) 4(fc+1)2 (fc-l) 4(fc - 1)2 



In 



(fc + 2) (fc+2) k^ 2 (fc-2) (fc - 2) 



for k < -2, 

for k = — 1 , 
for k = 0, 
for fc = 1, 

for fe > 2. 



In literature didn't find exact values of approximating coefficients. 
When q = 1 the Riesz-Feller operator is singular, hence the 
problem. Numerous works of Gorenflo and Mainardi i.e. |6*||7| 
propose various ways which determine values of the coefficients 
luj^*' (i.e. based on the Griinwald-Letnikov discretization) but 
they don't provide continuity in the interval a € (1,2]. The 
coefficients (23) can approximate the Cauchy process when we 
use (23) in numerical calculations. 

3.2. Fractional FDM 

While discretization of the Riesz-Feller derivative in space is 
done, in this subsection we describe the finite difference method 
for the equation of anomalous diffusion (1). Here we restrict this 
solution to one dimensional space. In comparison with the stan- 
dard diffusion equation where discretization of the second deriva- 
tive in space can be approximated by the central difference of sec- 
ond order, we will use generalized scheme given by formula (20). 
The differences appear in setting of boundary conditions. 

We shall introduce a temporal grid = t° < t 1 < . . . < 
t f < t /+1 < ... <oo with the step At = t /+1 - t f and we de- 
note value of the function C (x, t) in the point Xk at the moment 
oftimet / asC[ = C (x k ,t f ) for k eZand/ G N. 

3.2.1. Pure initial value problem 

In the explicit scheme of the FDM we replaced Eqn (1) by 
the following formula 



- C f 



At 

k — — oc 

After simplification finally we obtained 



oc) 



— E Ci+kPk 



where coefficients p£ are 

„.. .l + i^4 Q \f°rfc = 0, 
K a ^w[ a \ forfc/0. 



(24) 



(25) 



(26) 



Using simple calculations one may proof, that arise the fol- 
lowing relationship 



(27) 



In order to determine stability of the explicit scheme the co- 
efficient (26) for k = in formula (25) should be positive 



At 



,(«) 



(28) 



Hence we fixed the maximum length of the step At as 



At < 



fe°T(4-a) 



KaW M K a (23- - 4) (c L (a, 9) + c R (a, 9)) ' 



(29) 



The initial condition (9) is introduced directly to every grid 
nodes at the first step t = t°. This determines values of the 
function C as 



C° = c (xi) . 



(30) 



In unbounded domains the implicit method isn't easily ap- 
plicable because it generates infinite dimensions of all matrices. 
Thus one usually seeks improved difference equations within the 
explicit scheme. 

3.2.2. Boundary-initial value problem 

Presenting numerical solution (25) with included unbounded 
domain — oo < x < oo has no practical implementations in 
computer simulations. 

Now, we present solution of this problem on the finite domain 
£1 : L < x < R with boundary conditions (8). We divide this 
domain Q. into N sub-domains with h — (R— L)/N. Figure 1 
shows modified spatial grid. 

L R 
1 1 1 1 1 1 1 1 1 1 1 — | 1 — I 

— X.\ Xq X, Xn -^a'+i — 

Figure 1: The nodes grid over space 

Here we can observe additional 'virtual' points in the grid placed 
outside of the domain fl In order to introduce the Dirichlet 
boundary conditions we proposed treatment which bases on as- 
sumption that values of the function C in outside points are iden- 
tical as values in the boundary nodes xq or xn 



C(x hl t) 



G{x ,t)=g L {t), forfc<0, 
C(x N ,t) = g R (t), for fc> iV. 



(31) 



On the base of previous considerations we modify expression (20) 
for discretization of the Riesz-Feller derivative. Thus we have 



]T c{x l+k ,t)w[ a) 



+ 9 l (t) s L [ a) +g R (t)s R ^] 



(32) 



for i = 1, . . . , N — 1, where 

-i— 1 



(a) 

S Li — / ^ k 



(a) 
ID = 



T (4 - a) 



(33) 



OO 



3i J 



(a) 
w k = 



+ (l~lf- a ] C L , 



[-(N-i + 2) 3 



k=N-i+l 

3-a 



+3(N- i + 



T (4 - a) 
3 (N - if- a + (N - i - 1) 3 ~ Q ] c R 



(34) 



Putting this expression to Eqn (1) we obtain a finite difference 
scheme depending on weighting factor o. Here we assumed 



<?« ( At ( / + - 



(35) 
(36) 
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in order to simplify the numerical scheme. For internal nodes Xi, 
i — 1, . . . , N — 1 we have 



At 



= K a — 

h a 



N-i 



+9l s l 



k = — i 

(«) 



E ^/ +fc + (i-<7)c/+ 3 



/+! (a) 



and for the boundary nodes xq and a; at: 



C, 



c 



/+1 



with pj°"' defined by formula (26). 

We can observe that boundary conditions influence to all val- 
ues of the function in every node. In opposite to the second 
derivative over space which is approximated locally, the char- 
acteristic feature of Riesz-Feller and other fractional derivatives 
is dependence on values of all domain points. For a = 2 and 
(37) (9 = our scheme is identically as wide known and used the for- 
ward difference in time and central difference in space scheme 

(ftcs) rnigiin). 

The skewness parameter 9 has great significance influence on 



(38) 
(39) 



The method is explicit for a = 1 and partially implicit for 
< a < 1 and with a = being fully implicit. In literature this 
method is known as the cr-method for parabolic equations. 

Above scheme described by expressions (37)-(39) can be writ- 
ten in matrix form as 



A C f+1 = B, 
where 



(40) 



1 0. 

a_i 1 + ao ai a2 
a_2 a_i 1 + ao a± 
a_3 a_2 a-i 1 + a . 
a_4 a_3 a_2 a_i . 



a-JV+2 a-N+3 Cl-N+4 a-M+5 ■ 

a-N+i a-N+2 a-iv+3 a~N+4 ■ 
0. 





ajv-3 ajv-2 ajv-i 

ajv-4 ajv-3 ajv-2 

ajv_5 ajv_4 ajv_2 

ajv_6 ajv_3 ajv_4 

1 + ao ai a2 

a_i 1 + ao ai 

1 



(41) 



B 



f+i 
9l 

bi 

t>2 

b. 
In 



(42) 



OAT-2 

6jv_i 
f+i 

9r 

with 

a 3 = ( CT ~ 1 ) K "^ w j°'\ l ' 01 './ 

At 

+ ° E c Uk-! 



-N + 1,..., N -1,(43) 



6, = Cf + tf a 



/+5 (a) . /+§ (a) 

9l sl) +g R 2 s r ) 



if ( a ) 



k=—j 



for j = 1, 



(44) 



and C f+1 is the vector of unknown function's values C at the 
time t f+1 . 

Particular case of above scheme (37) is the explicit scheme 
(for cr = l) which may be simplified to 



Cl +1 = i 



f+h 

9 L : 



for i = 0, 



At ( f+l ( a ) /+! f tt ) 

[ 9l Ll 9r Rt 



9r 2 - 



+ E C4 fcP i Q) ,fori = l,...,iV-l, 

k — — i 



(45) 



for i = N, 



±1 + one can obtain the clas- 



the solution. For a — > 1 + and 
sical hyperbolic equation, i.e. the first order wave equation (the 
transport equation). In this case our scheme tends to the known 
Euler's forward time and central space (FTCS) approximation of 
Eqn (1). Unfortunately this is unconditionally unstable and there- 
fore this is disadvantage this method. 

Proposed numerical scheme makes a bridge between Gaus- 
sian and Cauchy processes. Our scheme is also a bridge between 
diffusion and transport phenomena. 

4. Simulation results 

In this section we present results of calculation. In all pre- 
sented simulations we assumed k a = lm a /s and the length of 
ID domain I = Ira. Figure 1 shows two charts over space (one in 
the logarithmic scale) with absorbing boundary C (x,t)\ x=0 = 
C (x, t)\ = 0. On these plots solutions for different values of 
parameter a G ( 1.01, 1.5, 2} at time t = 0, 0.01, 0.3s for = 
are presented. Figure 2 presents another example of the solution 



8 - 

u(x,t) 



/I 



T \ 
r f*\ 1 



( = 0s initial state 

after ( = 0.01.1 after ( = 0.3s 

— * — cc = 2 — ^ — a = 2 

— • — 0=1.5 — — 0=1.5 

— ■ — 0=1.01 — » — 0=1.01 



X 




Figure 2: Solution over space for a 6 ( 1.01, 1.5, 2) 
a) normal scale; b) logarithmic scale. 
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which differs from example presented by the Fig. 1 (boundary 
conditions C (x, i)\ x=0 = C (x,i)\ x=1 = 100 and initial con- 
dition C (a;, t) | t=0 = ). In both cases we observe diffusion 
process arising in different way. The last example reflects case 
when the parameter of skewness is 9 = 0.5 and a = 1.4. Fig- 
ure 3 shows a diffusion transport process over space at different 
moments of time. 




0,0 0,2 0,4 x 0,6 0,8 1,0 

Figure 3: Solution over space for a = 1.01, 1.5, 2. 
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Figure 4: Solution over space for a = 1.4 and 9 = 0.5. 



5. Conclusions 

In summary, we proposed the fractional finite difference met- 
hod for fractional diffusion equation with the Riesz-Feller frac- 
tional derivative which is extension to the standard diffusion. We 
analysed a linear case of diffusion equation and in the future we 
will work on non-linear cases. We obtained the implicit and ex- 
plicit FDM schemes which generalise classical schemes of FDM 
for the diffusion equation. Moreover, for a — 2 our solution 
equals to the classical finite difference method. 

Analysing plots included in this work, we can see that in the 
case a < 2 (the Levy flight) diffusion is slower then the standard 
diffusion (Brownian motion) in the initial time. Nevertheless, 
when we analyse the probability density function we observe a 
long tail of distribution in the long time limit. In this way we can 
simulate same rare and extreme events which are characterised 
by arbitrary very large values of particle jumps. 

Analysing changes in the skewness parameter 9 we observed 
interesting behaviour in solution. For a — > 1 + and for 9 — > ±1 + 
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we obtained the first order wave equation. For 9 6 (0, 1) (with re- 
strictions to order a) we generate a class of non-symmetric prob- 
ability density functions. 
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